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PACS 75 . 10 .Hk - Ising model, magnetic ordering 
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Abstract - We study the dynamics of 2d spin-ice following a quench from a fully disordered 
initial condition (equilibrium at infinite temperature) into its disordered, ferromagnetic and anti- 
ferromagnetic phases. We analyze the evolution of the density of topological defects and we show 
that these take finite density over very long periods of time in all kinds of quenches. We identify 
the leading mechanisms for the growth of domains in the ordered phases and we evaluate the 
(anisotropically) growing lengths involved in dynamic scaling. 



•^^ In a large class of condensed-matter systems the ten- 
^dency to local ordering is hampered by constraints. Frus- 

^tration entails, typically, a non- vanishing entropy at zero 
r^ temperature. 

I The prototypical example is water ice for which this zero 

^^ point entropy has been measured in the 30s 1 . Pauling 

^ explained this feature with a model in which the O atoms 

Q occupy the vertices of a coordination four lattice and ex- 

""Hactly two H atoms are near while the other two H atoms 

p^ are shifted away from each vertex J2] . The large degener- 

^ acy of such locally electro- neutral ground states gives rise 

00 to the zero point entropy. 

1^ A residual entropy has also been measured in frustrated 
^N^l magnets 3,4 such as Ho2Ti207 [5]. In these spin-ice sam- 
l^ples, the relevant interacting degrees of freedom are clas- 
^^ sical Ising spins located at the nodes of a pyrochlore lat- 
^H tice and aligned with the local axis connecting two corner- 
T-H sharing tetrahedra 4 . Each tetrahedron can be seen as 
^ a vertex in a Zd lattice taking one out of sixteen possible 
j configurations. The energy is locally minimized when two 
\ spins point inward and two outward verifying the ice rule 
^or zero-divergence condition for the magnetic moments 
carried by the spins. The other ten 'defects' have a mag- 
netic charge q = V - S. Experimentally, the statistical 
weights of the vertices can be tuned by applying uniax- 
ial pressure that break the isotropy of the ice-model but 
preserve the Z2 symmetry 6 . Magnetic fields along dif- 
ferent crystallographic axes can also be used to modify the 
vertex weights. 

Two dimensional spin-ice samples have also been pro- 






duced in the laboratory. On the one hand, high mag- 
netic fields project the system onto 2d Kagome slices 7 ; 
on the other hand, 'artificial' spin-ice on a 2d square lat- 
tice can be manufactured with litography [sIIq] . Although 
the magnetic moments in artificial spin-ice are a priori 
athermal, an effective thermodynamics can be generated 
with applied rotating magnetic fields [lO[[TTj. More re- 
cently, thermal excitations have also been engineered in 
these samples p!2l. 

In recent years, research in this field has been boosted by 
the exciting proposal that topological defects, in the form 
of magnetic monopoles and their attached Dirac strings, 
could be observed in spin-ice samples 13-18 . Following 



the proposal in |13] experiments have shown evidence for 
this kind of excitation 16-18 . Reaction-diffusion argu- 



ments have been used to estimate the density of defects in 
the disordered phase of 3d spin ice 15 . As far as we know 



no studies of dynamics in the ordered phases nor beyond 
this simple modeling have been performed. 

Here we choose a different approach to address the dy- 
namics of spin-ice models. For the sake of simplicity we 
focus on thermal quenches in the 2d square lattice spin 
ice model built as a stochastic extension of the celebrated 
6 vertex model of statistical mechanics [l9]. We use a 
reject ion- free continuous-time Monte Carlo (MC) algo- 
rithm 20 , with local spin- flip updates and non-conserved 
order parameter, that allows for thermally- activated cre- 
ation of defects. The longest time reached with this 
method, once translated in terms of usual MC sweeps, is 
of the order of 10^^ MCs, a scale practically unreachable 
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with usual Metropolis algorithms. This allows us to iden- 
tify the equilibrium phase diagram and to analyze different 
dynamic regimes. Our dynamic results are manifold. We 
reproduce known facts of the dynamics of spin ice samples, 
as the existence of long-lived metastable states although 
with no need of long-range interactions. We prove that 
the dynamics after a quench into the FM and AF phases 
conforms to the domain-growth scaling picture 121 and 
we identify the anisotropic growing lengths. We derive a 
large number of new results that, we propose, should be 
realized experimentally [22 . Our findings should also be 
of interest to the integrable systems community. 
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Fig. 1: (Color online.) The sixteen vertex configurations on 
the 2d square lattice and their weights. The first six vertices 
verify the ice-rule. The next pair completes the eight- vertex 
model and have charge ^ = ±2. The last eight vertices have 
charge ^ = ±1. This color code is used in Fig. [5] 

The model. We consider an L x L square lattice with 
unit spacing and periodic boundary conditions and we 
vary the linear size L from 10 to 300. Each edge is oc- 
cupied by an arrow modeled as a binary variable 6* = ±1. 
We assign a Boltzmann weight ujk oc e~^^^ to each of 
2^ four-arrow vertex configurations. The 

16 



the k = 1, . . . , 

energy is H = XlfcLi ^k^k^ where Uk is the number of ver- 
tices of type k. We set uji = UJ2 = a, ojs = UJ4 = 6, 
oj^ = ujQ = c for the ice-rule vertices and UJ7 = ujs = d^ 
ujg = ... = cjie = d for the 2-fold and 1-fold defects, re- 
spectively, ensuring invariance under reversal of all arrows 
(see Fig. [1]). Henceforth we measure the weights in units 
of c. 

Equilibrium properties. In the 2d six-vertex model de- 
fects are forbidden. A host of exact equilibrium properties 
of this model were derived with the Bethe Ansatz and via 
mappings to random matrix theory, algebraic combina- 
torics (domino tilings) and crystal growth 1 19 . Depending 
on the weight of the vertices the system sets into a quasi 
long-range ordered paramagnetic (FM) or disordered (D), 
two ferromagnetic (FM) and one anti-ferromagnetic (AF) 
phases separated by different transition lines. The four 
phases of the six- vertex model {d = 0) are characterized 
by the parameter Ae = (a^ + 6^ - l)/{2ab) [To]. (i) For 
Ae > 1 and a > b -\- 1 the vertices of type 1 and 2 (type 3 
and 4 for 6 > a + 1) are statistically favored and the sys- 
tem is frozen into one of its two symmetric ground states 
with perfect FM order and no low energy excitations. At 



|af ^ — b\ = 1, i.e. Ag = 1, there are two first-order phase 
transition, (ii) For — 1 < Ae < 1 the system is quasi long- 
range ordered. The ice-rule is strong enough to prevent 
full de-correlation at any temperature. At a^^ = 1 — 6, 
i.e. Ae = — 1 there is a Kosterlitz-Thouless (KT) phase 
transition, (iii) For Ae < — 1 vertices of type 5 and 6 are 
favored and the system orders into an AF state populated 
by low energy excitations. The phase diagram is shown in 
Fig. [2] with solid (red) lines. This scenario is modified by 
the ice-rule breaking vertices {d 7^ 0) as we discuss below. 
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Fig. 2: (Color online.) (a, 6)-plane projection of the sixteen- 
vertex model phase diagram. The solid (red) lines are the 
six- vertex model first-order transition between disordered (D) 
and ferromagnetic (FM) phases, and KT transition between 
disordered (D) and antiferromagnetic (AF) phase. The second- 
order transition lines for d / are shown with dashed (blue 
and black) lines. 

The 2d sixteen-vertex model is isomorphic to the Ising 
model on the checkerboard lattice with many-body inter- 
actions [23]. While this model is integrable for a special 
set of parameters for which the equivalent Ising model has 
two-body interactions only, none of these special cases cor- 
responds to our choice of parameters. We then used MC 
simulations to elucidate the equilibrium phase diagram 
and statistical properties of relatively small samples. We 
set the origin of coordinates on a vertex. The indices (a, /3) 
are the coordinates of each vertex and {{2a + l)/2),/3) 
and (a, (2/3 + l)/2) locate the mid-points of the right- and 
up- pointing bonds. We checked equilibration with stan- 
dard criteria; in particular, we verified (i) that the den- 
sity of vertices stabilizes at long times at the same value 
starting the evolution from very different initial configu- 
rations; (ii) that the spatially local two-time correlation 
function depends on the time difference only. These tests 
will be shown in [24]. After verifying equilibration in this 
way we computed: (i) The absolute staggered magneti- 
zation per spin defined as M± = ((|m^|) + (|m^|))/2 



with L'^ttVa 

and its counterpart m^ 



zl{o 



Sf^ 



■ Sf^ 



l3)eA ^(2a+l)/2,^±E(a,^)GB ^(2a+l)/2,/3 

^. In the expressions for m we di- 
vided the lattice into two sub-lattices A and B such that 
a -\- P is even and odd, respectively. The ± signs al- 
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Fig. 3: (Color online.) Study of the FM transition, (a) Magne- 
tization per spin M+ and (b) magnetization cumulant Km^ as 
a function of a for 6 = 0.5, d— 10~^ (the group of curves on the 
left) and d — 0.1 (the ones on the right), and several L given in 
the key. The crossing points of Km^ determine ac{b,d). The 
vertical dotted (black) lines are the critical values predicted by 
|Ai6| = 1. The horizontal dotted level is 1/3. 



low one to distinguish between FM and AF order, (ii) 
The fourth-order cumulant Km± = (^m^ +^m^ )/2 with 
Km^ = 1 - ((m^)^)/3((m^)^). (• • •) denotes an aver- 
age over independent runs. In Fig. p^ (a) we show M+ 
for b = 0.5 and two values of d as a function of a. 
L = 10, . . . , 40 and we averaged over 10^ — 10^ samples. 
The variation of the magnetization as a function of a shows 
a sharp jump at a = 1.5 for b = 0.5 and (i = 0, as expected 
for a first order phase transition. For d > the curve takes 
a sigmoid form that gets wider (less step-like) for increas- 
ing d. The intersection point appears at larger values of 
a for larger values of d. These features suggest that the 
transition to the FM phase is continuous, occurs at larger 
values of a, and that there are fluctuations in the ordered 
state. The crossing of Km+ at height 1/3 (dotted hor- 
izontal line) for several values of L shown in Fig. [3] (b) 
determines ac{b,d). Finite size scaling will be discussed 
in [24]. Consistently with a second-order phase transition, 
Km+ remains positive for all L and the energy cumulant 
(not shown) converges to zero. The analysis of the stag- 
gered magnetization and its cumulant demonstrates that 
the AF transition no longer belongs to the KT universality 
class and also becomes second-order as soon as d 7^ 24 . 
The projection of the critical lines onto the (a, 6)-plane 
are shown in Fig. ^ with dashed (blue and black) lines. 
For d = 10~^ they agree, within our numerical accuracy, 
with the exact critical lines for the six-vertex model. For 
increasing d the extent of the PM phase increases. The 



excitation properties are radically modifled by the relax- 
ation of the ice-rule: the paramagnetic (PM) phase loses 
its quasi long-range order and the FM state admits excita- 
tions |24 . These conclusions are in agreement with exact 
computations on the eight- [19, 23 and sixteen- 23 vertex 
models for special values of the parameters. We conjecture 
that the phase diagram is characterized by the anisotropy 
parameter Aie = [a^ + 6^ - c^ - {4d)'^]/[2{ab + c{4d))]. 
The PM phase corresponds to the parameter space sub- 
manifold with lAiel < 1, the FM phase to Aie > 1 and the 
AF phase to Aie < —1 in agreement with the fact that the 
transition lines are parallel to the six- vertex model ones. 
A similar displacement of the critical lines was found in 
spin-ice on the Husimi tree 125 . 




Fig. 4: (Color online.) Time-dependent density of defects, 
nd{t), after a quench from T^ootoa = b=l and d = 10~^, 
10~^ . . . , 10~^. (a) L = 50 and (b) L = 100. The black curves 
are ford^lO"^, 10~^, 10~^. The gray (color) curves are for 
smaller values of d decreasing from left to right, (c) Short time 
behavior in the case d = 10~^ and L = 100 confronted to the 
decay po/(l + ^t) (dashed curve) [15 and the fit po/{l-\- Qt)"" 
with a = 0.78 (blue plain curve), (d) Test of scaling with td^ 
for systems with L = 50. 

Quench dynamics. We now turn to the dissipative 
stochastic dynamics of an equilibrium initial conflgura- 
tion dita = b = d=l (i.e. T -^ 00) after a quench to 
sets of parameters in the (i) disordered, (ii) FM, and (iii) 
AF phases. In case (i) the system should equilibrate easily 
but the question remains as to whether it gets blocked in 
metastable states with a large density of defects. In cases 
(ii) and (iii) the interactions between the spins, mediated 
by the choice of vertex weights, should create ordered do- 
mains, FM or AF. The quantitative characterization of 
growth in the ordering processes is given by two possibly 
different growing lengths extracted from correlation func- 
tions along orthogonal directions || and ± that we identify. 
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Fig. 6: (Color online.) Interfaces between FM domains. Local 
vertices and spins on the bonds are shown. Left panel: diagonal 
wall (red solid line). Central panel: a 'loop' fluctuation on 
the plaquette highlighted in the left panel. Right panel: a b 
corner vertex cannot be neighbor of an a- vertex, explaining the 
presence of strings. 




Fig. 5: (Color online.) FM ordering. Upper panel: time evo- 
lution of the density of vertices with weight a,b,c,d for a = 5, 
b = 1, d = 10~^ and L = 100 averaged over 300 samples. The 
snapshots are typical configurations at the dates indicated by 
the arrows. Black/white points are vertices 1/2 and the rest 
are shown in gray (color) scale. Inset: time-dependence of the 
longitudinal (upper curves) and transverse (lower curves) grow- 
ing lengths for three system sizes, L = 100, 200, 300. A fit to 
t'^^^ is shown with a dotted black line. 



(i) Quench into the disordered (D) phase. Figured dis- 
plays the time-dependent density of defects, nd{t)^ de- 
fined as the number of vertices of type 7-16 divided by 
I/^, after an infinitely rapid quench to a = 6 = 1 and 
d = 10~^, . . . , 10~^ of samples with linear size L = 50 (a) 
and L = 100 (b). These data have been averaged over 10^ 
runs. 

For large d (black dark curves) nd{t) quickly saturates 
to its equilibrium value. Numerical estimates of the equi- 
librium density of defects, n^^ for d = 10-\ 10-^ 10"^ 
are shown with dotted black lines. As expected n^f is 
an intensive quantity that increases with d. It does not 
depend upon the system size for L > 50 and d > 10~^. 

For small d (~ 10 ~^) the systems do not reach equilib- 
rium within the simulated time- window. After a first de- 
cay, nd{t) gets frozen at approximately constant values be- 
fore relaxing, in a much longer time-scale, to a configura- 
tion in which only two defects are present in our small sam- 
ples. Note that in order to distinguish the (i— dependent 
equilibrium values for these very small ds one would need 
to equilibrate much larger samples. Unfortunately, the 
special purpose loop algorithm devised for the 6 vertex 
model does not apply to our generalized case. For our 
working sizes we see asymptotes taking the value 2/L'^ on 
average. 

The initial decay is fitted by a po/(l + ^^)" power-law 



decay with a c^ 0.78 over three orders of magnitude in 
t and Udj as shown in panel (c) in Fig. [4] The power- 
law is shown with a solid blue line in the figure together 
with the data for d = 10~^ and L = 100. This law is 
different from the simple t~^ decay found with a mean- 
field approximation to a diffusion-reaction model shown 
with a dashed red line in the same figure 1151 . 

The decay is next arrested at a metastable density of de- 
fects n^ ^ 10/L^. The plateau lasts longer for smaller d 
and its height is roughly independent of d. This feature is 
reminiscent of what was found numerically in dipolar spin- 
ice although contrary to the modeling in fl^ our model 
does not have long-range interactions. At the entrance to 
the plateau the system has between 3 and 4 times more 
defects of type 7-8 than those of type 9-16 which are only 2 
or 3. Therefore, in the final decay from the plateau to the 
asymptotic value rid ^ 2/L^ the remaining doubly charged 
defects have to disappear. This may be due to two kinds 
of processes. In the first case, two defects of type 7 and 8 
meet to produce two singly (and oppositely) charged de- 
fects with no energetic gain. The total density of defects 
remains constant after this reaction. An example of the 
second case is a reaction in which a defect of type 7 (charge 
q = 2) meets one of type 14 (charge q = —1) to produce 
a defect of type 10 (charge g = 1) and a spin-ice vertex 
with no charge. In terms of a reaction-diffusion model this 
corresponds to 2q + {—q) ^ q -\- and, consequently, an 
energetic gain. Note that the number of single charged de- 
fects has not been modified in this process but the number 
of doubly charged defects diminished and so did the total 
number of defects. In both cases the remaining defects 
need to diffuse, a process with no energetic cost, to find 
a partner and annihilate. From inspection of the individ- 
ual runs and the densities of single and doubly charged 
defects we see that the second process is favored, as also 
suggested by the energetic gain. This regime is character- 
ized by a scaling of the dynamic curves with the scaling 
variable td^ 26 as shown in Fig. H (d) for the L = 50 
data. 

(ii) Quench into the FM phase. We choose a = 5, 6 = 1 
and d = 10~^, favoring vertices with weight a. In Fig.|5]we 
present the density of vertices, n^{t), with n = a^h^ c, d, in 
a log-linear scale. The evolution is illustrated with three 
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Fig. 7: (Color online.) Schematic representation of FM stripe 
motion. Vertices on each site are specified. Diagonal (red) 
lines delimit domains of opposite magnetization. Black arrows 
indicate the spins that flip to get the new configuration (rep- 
resented in blue after the flip). 



configurations at instants shown with vertical arrows. Do- 
mains grow anisotropically and we choose the || and ± di- 
rections to be parallel and perpendicular to the diagonal 
joining the lower-left and upper-right corners in the pic- 
tures, respectively. During a short transient {t ~ 0.01 
MCs) all densities remain roughly constant (regime I). 
Suddenly, a large number of defects are transformed into 
divergence- free vertices by a few single spin- flips: rid de- 
cays while n^, 77,5 and Uc increase (regime II) indepen- 
dently of a 24 . A typical configuration at this stage is 
the left-most snapshot and there is no visual ordering as 
corroborated by the small values taken by Ly ^ and dis- 
played in the inset in a log-linear scale for three values of 
the system size, L = 100, L = 200 and L = 300. Sub- 
sequently the system sets into a slow relaxation regime 
in which the dominant mechanism is the one of growing 
anisotropic domains with FM order, see the central snap- 
shot (regime III); n^ depend upon a and there are as many 
domains with m^^y = 1 (vertices 1) as m^^y = — 1 (vertices 
2) respecting symmetry. In this regime Ly grows faster 
than L± and tends to saturate to an L-dependent value 
when the stripes are fully formed. For the largest sample 
size, L = 300, our numerical data are consistent with a 
^1/2 growth that is shown with a dotted black line. In- 
stead order in the ± direction has not yet percolated. The 
full equilibration of the sample needs the percolation of 
order in the ± direction which is achieved by a still much 
slower mechanism (regime IV). 

A better understanding of the processes involved in the 
ordering dynamics is reached from the analysis of the snap- 
shots, (a) Domain walls are made of c-vertices and pla- 
quettes of divergence- free vertices, as shown in the left and 
central panels in Fig. |6) respectively. The latter are 'loop' 
fluctuations in which all the spins on the plaquette are 
sequentially flipped. Interfaces between FM states tend 
to be parallel to the main diagonal, which one depending 
on which FM phase one quenches into, (b) Quasi-one- 
dimensional paths made of b and c vertices (loop fluctua- 
tion can be attached to them) act as bridges between two 
domains of the same type and run through a region with 



the opposite order. These structures are similar to the 
ones found in the kinetically constrained spiral model 27 . 
In order to further increase the density of a- vertices and 
develop the FM order the domain walls and bridges have 
to be eliminated. The latter disappear first via the follow- 
ing mechanism. 'Corners' made of b (or, less commonly, 
d) vertices sit on a curved domain wall. Such b vertices 
cannot be surrounded by more than two type 1 or 2 ver- 
tices (only defects can, see the third panel in Fig.lGJ). The 
string progressively disappears eaten by the attached do- 
mains that grow from the corner or, alternatively, it is 
first cut by the creation of two defects and the two strands 
subsequently shrink, an extremely slow process. Once the 
path has been eliminated one is left with two defects sit- 
ting on the walls of the now detached domains, that move 
along the interface and eventually annihilate with their 
anti-partner. Once parallel bands are created (third con- 
figuration in Fig. Is]) the mechanism in Fig. [tI takes over 
(regime IV). After the creation of a pair of defects on the 
interface, the sequence of steps in the figure shrink the 
vertex 1 stripe on a time scale that diverges with L. 




Fig. 8: (Color online.) AF ordering, 
density of vertices in a system with L ■ 



Time evolution of the 
= 100 after a quench to 



a = 0.1, 6 = 0.1, d= 10~ averaged over 300 runs. Inset: the 
time-dependent growing length L\\ confronted to t^^^ (dotted 
black line). Typical configurations are shown. 



(iii) Quench into the AF phase. The evolution of the 
vertex population is shown in the main panel in Fig. |8] 
for a = 6 = 0.1 and d = 10~^. Similarly to what found 
in the FM quenches, in regime I all densities remain ap- 
proximately constant. This is followed by regime II with a 
rapid annihilation of defects into divergence-free vertices. 
The creation of a, b and c-vertices occurs with a rate that 
depends on a while, surprisingly, rid does not, at least 
within our numerical accuracy. In regime III the system 
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increases the AF order by growing domains of staggered 
magnetization ±1 with c vertices. Since a is very close to 
h for our choice of parameters, domains are quite isotropic 
and the growing length are, within numerical accuracy, 
t^/^. Regime IV follows next and it is characterized by 
a strong slowing-down although there is no obvious ex- 
tended structure blocking the evolution. In regime V the 
system finally reaches equilibrium. The relevant elemen- 
tary moves in each regime will be discussed in [24| . 

We presented a numerical study of the quench dynam- 
ics of 2d spin-ice. In this Letter we established the main 
statics and dynamical properties of this model: equilib- 
rium phase diagram and relaxation dynamics after ther- 
mal quenches. No study of frustrated magnet existed in 
the litter ature such complete. We demonstrated the ex- 
istence of long-lived metastable states with an excess of 
topological defects in a model with no long-range inter- 
actions, cfr. 15 . We showed that the dynamics after 
quenches into the FM and AF phases proceed by coarsen- 
ing of equilibrium domains. More details will be presented 

Our results should mo- 
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in an extended publication 
tivate new experimental studies of the dynamics in frus- 
trated magnets and possibly motivate researchers in the 
integrable systems community to try to adapt their very 
powerful techniques to deal with dynamical issues. 
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